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We have studied boron nitride monolayer and bilayer band structures. For bilayers, the ground 
state energies of the different five stackings are computed using DFT in order to determine the most 
stable configuration. Also, the interlayer distance for the five different types of stacking in which 
boron-nitride bilayers can be found is determined. Using a minimal tight binding model for the 
band structures of boron nitride bilayers, the hopping parameters and the onsite energies have been 
extracted by fitting a tight binding empirical model to the DFT results. 
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INTRODUCTION 



It was shown, since the early days of graphene re- 
search, that other layered materials can be exfoliated 
by micromechanical cleavage, following exactly the same 
procedure used to isolate graphene from graphite, giving 
birth, in this way, to new two-dimensional crystals, which 
will trigger, with certainty, new scientific investigations»i 
Among these new materials there are several dichalco- 
genides, complex oxides, and boron nitride. It was also 
shown that these new two-dimensional forms of con- 
densed matter show high crystal quality and macroscopic 
continuityji 

This new world of two-dimensional crystals is still 
much unexplored, both in what concerns the fundamental 
properties of these systems and their potential applica- 
tions. Chief among these two-dimensional crystals, and 
excluding graphene for obvious reasons, boron nitride is 
emerging as a system which is finding applications as a 
scaffold for graphene in electronic devices. The fact that 
boron nitride is an insulator, having a large energy gap 
(^5 — 7 eV), together with its high degree of purity makes 
it a perfect membrane to isolate graphene from the im- 
perfect gate dielectrics used so far. Since both systems, 
that is, graphene and boron nitride, interact weakly, two- 
dimensional crystals of the latter can be used as a buffer 
layer between a substrate and graphene, leaving its in- 
trinsic electronic properties unafected by disorder associ- 
ated with the dielectric surface. As a consequence, large 
electronic mobilities, of the same order of magnitude of 
those found in suspended graphene samples, have been 
measured in graphene on top of boron nitridei^ Graphene 
devices on hexagonal BN (ft,-BN) have been shown to 
exhibit enhanced mobility, reduced carrier inhomogene- 
ity, and reduced intrinsic doping, in comparison with 
graphene laying directly on top Si02'^ Additionally, this 
new experimental setup allowed the measurement of the 
fractional quantum Hall effect using the four probes ge- 
ometry, exposing a buried panoply of fractional filling 
factors, still waiting for deep fundamental justification. 
Other applications of BN can be envisioned, ranging from 
spacers to tunneling barriers. 



Since these first applications of BN in fundamental 
research that the interest in single layer boron nitride 
(sBN) has been increasing steadily. In the recent past 
several experimental studies have been performed)^^— 
including TEM^"— and AFM^ characterization of the 
material's surface, optical and Ramani spectroscopy, 
and intentional damage sBN's surface.^^ The main ex- 
perimental approaches for producing sBN have been 
both mechanicali'^'^i^ and chemical^ exfoliation. Oth- 
ers chemical methods are also available nowi^i^ In addi- 
tion to the discovery that sBN is an excellent substrate 
for graphene electronic devices^ the fabrication of large 
area (several cm^) BN layers?^ opens the possibility of 
mass production of this new two dimensional material, 
at the time of writing quite an expensive one. 

In parallel to the experimental investigations, some 
theoretical studies on sBN and bilayer boron nitride 
(bBN) have been performed; Topsakal et al^ made a 
density functional theory (DFT) study of the electronic, 
magnetic, and elastic properties of sBN and boron ni- 
tride nanoribbons. Giovannetti et alr^ studied graphene 
on the top of ft,-BN by DFT. These authors concluded 
that graphene opens a small gap for three orientations of 
graphene on the top of /i-BN, and that the most stable 
orientation was an AB stacking with the carbon atoms 
on the top of the boron atoms. Following the investi- 
gations by Giovannetti et al., Slawinska et alr^ studied 
graphene AB stacking on a sBN, using both tight-binding 
(TB) and DFT. They confirmed the opening of the gap 
and fitted the TB parameters to the DFT calculations. 
The difficulty with the aforementioned calculations is the 
fact that they assume graphene and sBN unit cells of the 
same size and perfectly oriented. We note in passing that 
the opening of a gap in graphene's spectrum, when this 
material is laying on top of sBN, has been speculated to 
exist in the past and the transport properties of graphene 
under these conditions have been computedJ^ 

Despite the above theoretical studies, it is however ex- 
perimentally known that graphene on the top of sBN has 
in general a random crystallographic orientation, that is, 
there is no preferential orientation of graphene's lattice 
relatively to that of sBN<^ In addition, no energy gap in 
graphene's spectrum has been measured so fari^ 
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FIG. 1. (Color online) Illustration of the honeycomb lattice 
with the A and B sublattices, the unit cell, and the primitive 
vectors ai and 02. It is equivalent to assume that boron or 
nitrogen is on the A or on the B sublattice. 



Marom et al— studied the interlayer sliding of one 
layer of BN on the top of another layer of BN using DFT 
with van der Walls corrections. They found that the 
lower energy configurations (AA' and AB) of bBN differ 
in energy by an amount smaller than the accuracy limits 
of their calculations. 

In this paper we investigate the structural nature of the 
ground state of boron nitride bilayers using density func- 
tional theory. We consider the different possible stack- 
ing of the individual BN planes, compute their energy 
ground states and band structures, and fit the TB hop- 
ping parameters and on-site energies to the DFT results. 
Graphene bilayer can have two types of stacking: AA and 
AB (Bernal stacking), with the AB structure being the 
most stable one. BN bilayers have richer structural possi- 
bilities; given that there are two types of atoms, there are 
five possible bilayer stacking: two AA and three AB. Ex- 
periments show that both types of stacking exist when 
one considers multilayer BN.— It is then important to 
clarify what can be expected in what concerns the ener- 
getic stability of the different stacking when one takes a 
bilayer as an isolated system. 



II. A MINIMAL TIGHT-BINDING MODEL FOR 
BILAYER BORON NITRIDE 

Since we have in mind possible applications where both 
graphene and boron nitride two-dimensional crystals are 
brought into close contact with each other, we will be 
most interested in the energy bands due to the hybridiza- 
tion of the Pz orbitals alone, the so called tt— bands. The 
TT— bands of BN can then be studied using an empirical 
tight-binding approach. The goal is simple: using DFT 
methods we will parametrize the tight-binding hopping 
and onsite energies of boron nitride layers, which can lat- 
ter be used for microscopic calculations, such as tunneling 
properties. 

In Fig. [T] we represent the unit cell of a single layer 
of boron nitride with the primitive vectors ai and 02, 
both of length a. The minimal tight binding model for 
a sBN has three parameters only: the hopping t among 
nearest neighbor atoms and the onsite energies at the 
boron and the nitrogen atoms (in fact, one of these two 



onsite energies can be taken as a reference energy, and 
one ends up with two fitting parameters only) . According 
to this simple model, the Hamiltonian for the electrons 
in the tt— bands of sBN can be written as 



where Eb is the energy at the boron site. En is the 
energy at the nitrogen site, and 

(j)/t=l + gia(-'=x/2+V3fc„/2) _|_ gM(fc,/2-|-V3fc„/2) _ ^2) 

The wave vector in the Brillouin zone is written as fc = 
(kx, ky) and the eigenvalues of Hamiltonian ([T|) are given 
by 

E^Eo±]^^El+m\ (3) 

where Eq = {Eb + En)/2 is the energy in the middle of 
the gap and Eg = Eb — E^ is the energy gapjii 

Using the above defined tight binding model one has 
three parameters that one needs to adjust: the hopping 
parameter t, the band gap Eg, and the middle gap energy 
-Bo, which can conveniently be chosen to be zero. These 
parameters can be found by adjusting Eq. ([3]) to the 
band structure of sBN computed from first principles. 

For bilayer boron nitride, our minimal model includes, 
in addition to the one given above, an interlayer hopping 
parameter t' between the two atoms that arc directly 
on top of each other. Hopping between atoms located 
at larger distances can be added to the model, at the 
expenses of having more fitting parameters. If one wants 
to fit accurately at the same time the valence and the 
conduction bands, adding these extra hoppings is strictly 
necessary. However, it is known that conduction bands 
in insulators are not well described by DFT calculations 
and thus we have found not to be necessary to go beyond 
this minimal model. 

Five different structures of BN bilayers can exist in 
principle, assuming that either the atoms are located ex- 
actly on top of each other or at the center of the hexagons 
(we are excluding, for example, twisted bilayer&i^); the 
five structures are represented in Fig. [2] (in the text that 
follows we will be referring to the different structures of 
BN bilayers by the notation introduced in that figure, 
which refer to the possible types of stacking). The AA 
and AA' stacking (see Fig. [2]) differ on what type of 
atoms are superimposed: in the AA case, atoms of the 
same type are superimposed, while in the AA' case the 
boron atoms arc on top of the nitrogen atoms and vice- 
versa. The A'B stacking has the nitrogen atoms superim- 
posed on the two layers and a boron atom at the center 
of the hexagon, while AB' stacking has the boron atoms 
superimposed on the two layers and a nitrogen atom at 
the center of the hexagon. The AB stacking has atoms 
of different types superimposed on the two layers and a 
boron atom at the center of the hexagon in one of the lay- 
ers, while on the other layer is a nitrogen atom that is at 
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FIG. 2. (Color online) Drawings of the five possible BN bi- 
layers. Some perspective was allowed to better recognize the 
type of bilayer. 



the center of the hexagon. The differences among these 
fives possibilities should be clear from a close inspection 

of Fig. m 

For boron nitride bilayers one finds that only two differ- 
ent Hamiltonians need to be written, since they describe 
all the five possible stacking (see Fig. [2|) of bBN, after 
the appropriate changes. Therefore we have: 



Haa 



El (j) t' 

<P* E2 t' 

t' E3 <j) 

t' (j)* E4 



(4) 



for the first two stacking in Fig. [2](AA and AA'), and 

(5) 



H 



AB 



El (j) t' 

cj)* E2 

t' E3 (I) 

(/)* E^ 



for the remaining three cases in the same figure (AB, 
AB' and A'B). In the last two Hamiltonians Ei, E2, E3 
and E4 are the energies at sites 1 to 4, two of which 
are boron atoms and the other two are nitrogen atoms. 
Atoms in the same plane are of different type implying 
that El ^ E2 and E3 ^ E4. The parameter t' is the 
hopping between planes and is the same expression as 
in Eq. [2] Each of these Hamiltonians give four bands. 

For the Hamiltonian (|4|), cases AA and AA', the four 
energy bands can be written as: 



E 



1 



Eo±yE^^+A{t'^ + m±8m, (6) 



where Eq and Eg have the same definition as above. For 
the Hamiltonian ([5]) we have to consider two different 
cases. The first one when Ei = E3 and E2 = E4, mean- 
ing that there are two atoms of the same type (belonging 



to different layers; cases AB' and A'B) on top of each 
other; in this case the energy eigenvalues read: 



t' 



1 



E = Eo±- + -^{Eg±t')^ + A\<l>\\ 

2 2 (7) 



t' 



1 



E = EQ±---^{Eg±t'f 



Note that Eg = Eb — En for the case AB' and Eg = 
En — Eb for the case A'B. The second case happens 
when El = E4 = Eb and E2 = E^ = E^, where the 
atoms on top of each other are of different types; in this 
case the energy eigenvalues are given by: 



1 



E^EQ±-^El±2t'^ 



-4|(/)|2 + 2i'2 + 4|0|2. (8) 



In equations ([6]), ([7]), and ([8]) there arc four parameters 
which one needs to adjust: Eq, which can be set to zero, 
Eg, t, and t' . In the bilayer. Eg is not the energy gap, 
although it has the same definition as in the monolayer 
(see above). Having described the TB model for bBN, 
it is necessary to perform first principle calculations in 
order to determine the bands of these systems. Once 
this is done the TB bands are fit to the ab-initio results 
and the TB parameters defined above are extracted. 

The density functional theory calculations were per- 
formed with an ab-initio spin-density functional code 
(aimpro)J^ We have used the GGA in the scheme 
of Perdew, Burke, and ErnzerhofJ^ Lower states (core 
states) were accounted for by using the dual-space sep- 
arable pseudopotentials by Hartwigsen, Goedecker and 
Hutterfi^ The valence states are expanded over a set of 
s— , p—, and d— like localized Cartesian-Gaussian Bloch 
atom-centered functions. In the framework of DFT calcu- 
lations, the Brillouin-zone (BZ) was sampled for integra- 
tions according to the scheme proposed by Monkhorst- 
Packjii The fe-point sampling was 20 x 20 x 1 and both 
the atoms and the unit cell parameter were relaxed iter- 
atively in order to find the equilibrium positions and size 
of the primitive cell. A supercell with hexagonal sym- 
metry was used; the parameter a was varied to find the 
equilibrium lattice constant, while the c parameter was 
kept at 60 a.u. 



III. DFT RESULTS AND TIGHT BINDING FITS 

OF THE ENERGY SPECTRUM OF 
MONOLAYER AND BILAYER BORON NITRIDE 

Boron nitride monolayers have very different electronic 
properties from those of graphene, since the broken sym- 
metry of the A and B sub-lattices necessarily precluded 
the existence of Dirac cones at the corners of the Bril- 
louin zone, creating a large band gap (larger than 5 
eV) in the BN band structurci^ Except for this vital 
difference, the rest of the electronic band structure of 
sBN resembles that found for graphene. The gap is at 
the isT-point and not at the F-point (as in usual semi- 
conductors) and, except for that particular zone in k- 
spacc, the shape of the bands of BN is similar to that of 
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FIG. 3. (Color online) Comparison between the band struc- 
ture of a BN monolayer and that of graphene. The solid lines 
refer to the bands of BN and the dotted curve to the bands 
of graphene. As discussed in the text, the main difference 
between the two band structures occurs close to the K and 
K' points. 



graphene bands. This is clearly shown in Fig. [3l where 
the band structure of the boron nitride single layer is 
shown superimposed on the band structure of graphene 
for comparison. The lattice parameter we have obtained 
for BN \s aBN — 0, = 2.51 A, in good agreement with 
the value measured experimentally?^ and for graphene we 
have found Qg = 2.46 A, giving a lattice mismatch of the 
order of 2% between the lengths of the primitive vectors 
of both crystal structures. 

The parameters of Eq. ([3]) were adjusted to the tt 
orbitals of the valence band only, since DFT is more ac- 
curate for occupied states. The fitting was done to all 
the band \tl th.e K - T - M - K' line, unlike the fit 
done in Ref. [ll|. In this latter work the bands were 
fit near the -point only. The curves resulting from 
fitting Eq. ([3]) to the ab initio data are shown in Fig. 
m The fitting to the valence band is excellent, but the 
resulting parameters do not give a perfect fit to the con- 
duction band, although the TB conduction band does 
follow the same trend as that obtained from the ab-initio 
calculation, as it should be. Given that DFT calcula- 
tions results in conduction and valence tt— bands that 
are not particle-hole symmetric (for example, one has 
the width of the conduction band greater than that of 
the valence band), and that Eq. ([3]), by construction, 
preserves that symmetry, the tight binding model con- 
sidering only a nearest neighbor hopping parameter will 
never fit both bands at the same time. The tight-binding 
parameters obtained from fitting the ab-initio bands are 
Eg = 3.92 eV and t = 2.33 cV, and the best fit is plot- 
ted in Fig. (|4]). The number we have obtained for t is 
different from the one (t=2.79 eV) obtained in Ref. [ll| 
and also from that obtained by the same group in a more 
recent work (t=2.28 eV) In these two works the fitting 
was done only near de iiT-point, and considering both 
valence and conduction bands, which is, most likely, the 
reason for the discrepancy between ours and the latter 
resultiS by those authors. We note that in their most re- 
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FIG. 4. (Color online) DFT band structure of the tt— bands 
of a BN single layer (solid line) and of the tight binding fit 
(dotted line). As explained in the text, only the valence band 
was fit. 



TABLE I. Interlayer distances for the bilayer samples, as cal- 
culated by DFT. The central column specifies the two atoms, 
each of a different layer, that are considered for the measured 
distance. 



Sample 


Atoms 


d(A) 


AA 


B-B 


3.75 


AA 


N-N 


3.75 


AA' 


B-N 


3.57 


AB 


B-N 


3.57 


A'B 


N-N 


3.72 


AB' 


B-B 


3.60 



cent worki^ those authors have found a value for t closer 
to the one we are reporting here. Our result for t also 
agrees very well with the one we have obtained for the 
bBN (see Table IIIip : this is an important consistency 
check for our calculations. 

We now turn to the study of the electronic spectrum 
of boron nitride bilayers. In our DFT calculations the 
primitive vector length a was let to relax for all stack- 
ings and its final value was found to be a = 2.51 A for 
all the cases represented in Fig. [21 the area of primitive 
cell size remains essentially unchanged for the five types 
of stacking considered here. The situation is different in 
what concerns interlayer distances, as expected. Table |l] 
summarizes the results: if different types of atoms are su- 
perimposed, the two planes are closer than if the atoms 
are of the same type, since the chemical bonds of the 
type boron-nitrogen are preferred (leading to lower en- 
ergy states), whereas boron-boron and nitrogen- nitrogen 
bonds are not. The two most stable stackings (A A' and 
AB) have the same interplane distance. The interplane 
distances calculated here are consistent, although some- 
what larger, than the experimentally measured value, 
which is - 3.3 A.i^i^ 

Taking the most stable stacking (lowest ground state 
energy) as a reference (AB stacking), the energy differ- 
ence (per unit cell) of the corresponding ground states 
of the different stackings are given in Table |lTl The 
AA stacking is the most energetic, followed by the A'B 
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TABLE II. Energy difference, per unit cell, relatively to the 
most stable stacking AB and the band gaps Ag for the five 
bilayer systems, as calculated by DFT. The third column in- 
dicates whether there is a degeneracy of the 7r-bands at the 
K and K' points in the Brillouin zone (see Fig. [5]). In the 
last column c.b. and v.b. stand for conduction and valence 
bands, respectively (see also Fig. [SJ. 





E (meV) 


A, (eV) 


degeneracy 


AA 


13.5 


4.23 


yes 


AA' 


0.4 


4.69 


no 


AB 


0.0 


4.60 


no 


A'B 


10.5 


4.52 


yes (c.b.) 


AB' 


3.0 


4.29 


yes (v.b.) 



AA' 



• GGA 
— TB 



K r M K' 



AB AB' A'B 




K K K 



FIG. 5. (Color online) Electronic band structure of boron ni- 
tride AA' bilayer and details of the electronic band structure 
near the K-point for BN bilayers AB, AB' and A'B. The solid 
line stands for the tight binding fits and the dotted lines for 
the ab-imto band structures. We note that only the valence 
bands have been fit, for the reasons given in the text. We can- 
not expect therefore that the conduction bands are accurately 
described by the tight binding energy bands. 



one, in which case the nitrogen atoms are on the top of 
each other. It is interesting that the most stable stack- 
ing is not the AA' stacking as found for bulk BNf2£, but 
the AB stacking instead, as found for bilayer graphene 
and graphite. These results are consistent with those 
obtained in Ref. [l^. As in that reference, the differ- 
ence between the ground state energies of the AA and 
A'B stacking is very small, meaning, in practice, that 
both types of stacking can exist; in fact, both have been 
observed cxperimentally<^ Table |ll] also shows the band 
gaps for all the bilayers band structure. The differences 
among them are small, reflecting the similarity of the 
band structure of these systems, both qualitatively and 
quantitatively. The differences in the band structure of 
the different types of bilayers are only noticeable near 
the K and K' points, as can be seen in Fig. [5] Figure [5] 
shows the full band structure of AA' bilayer and the de- 
tails of the electronic band structure near the -point for 
the others bilayer stackings, ploting both the DFT cal- 



TABLE III. Tight-binding parameters for the different bilayer 
stacking. The negative Eg value for A'B results from the 
exchange between N and B atoms for that particular stacking. 
In the table. Eg, t and t' stand for the Eb — En {En — Eb 
for A'B), the intralayer, and interlayer hopping parameters, 
respectively. 





Eg (CV) 


t (cV) 


t' (eV) 


AA' 


4.08 


2.36 


0.32 


AB 


4.19 


2.37 


0.60 


A'B 


-3.91 


2.34 


0.25 


AB' 


4.32 


2.38 


0.91 



FIG. 6. (Color online) Isosurface for the electron density 
in the 7r-band for the BN bilayer A'B. It can be seen that 
the electronic density is much higher at the nitrogen atom 
(smaller spheres) than at the boron atoms (larger spheres). 

culations and the tight binding fits. The tight-binding 
parameters for each bilayer represent in Fig. ^ Eg, t and 
t', were adjusted such that Eqs. ([6|), jT]), and ^ gave 
the best possible fit to the corresponding band structures 
as calculated by DFT. Only the valence band fits are ac- 
curate, since they are more reliable on DFT calculations 
than the conduction bands (a consequence of the fact 
that conduction bands have no electronic density in in- 
sulators). The TB conduction bands, as resulting from 
the fitting to the valence bands, are not in quantitatively 
good agreement with the DFT calculations, but qualita- 
tively do show the same features. The resulting param- 
eters arc shown in table IIIII It can be seen in Fig. [5] 
that the fits to the ab-initio bands are well described by 
the tight binding model using the parameters given in 
Table Hill As in the case of boron nitride monolayer, the 
tight binding valence bands describe quite well the bands 
calculated by DFT. 

Although the three bottom panels of Fig. [5] may not 
give a clear impression of the good quality of the fit of the 
valence bands, given the scales of the figures, the fits are 
as good as that seen in the top panel of the same figure, 
where the full band structure of the AA' bilayer is shown. 
For this reason, the full band structure is plotted only 
for the AA' stacking, while for the other cases only the 
details of the band structure close to the K' are shown. 

The values found for the intralayer t hopping param- 
eter are very similar for the different stackings, as are 
the onsite energies at the boron and nitrogen atoms; this 
is expected since the intralayer hopping should be about 
the same for all stackings as well as for a BN monolayer, 
and indeed it is. 
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The TT valence and conduction bands are degenerated 
at the ii'-point for the AA' bilayer. For the AB stack- 
ing, there is a shght energy difference (at the K-pomt) 
between the two tt— valence and conduction bands, and 
therefore no degeneracy is observed. For the A'B bi- 
layer the TT— conduction bands are degenerated at the 

-point, whereas for the AB' bilayer the degeneracy oc- 
curs for the TT— valence bands. It can then be concluded 
that when a boron atom is on the top of another of the 
same species, the tt— valence bands are degenerate at the 
JC-point, whereas when a nitrogen atom is on the top 
of another of the same kind, the reverse happens and it 
is the TT— conduction bands that are degenerate. It then 
comes with no surprise that for the AA' bilayer both the 
conduction and valence bands are degenerate, while for 
the AA bilayer they are not (these bands are not shown 
since they refer to the bilayer with (by far) the highest 
ground state energy when compared to the most stable 
bilayer). 

The electron density is higher around the nitrogen 
atoms than near the boron atoms (see Fig. [S]), which 
results in a larger distance between the two planes of the 



bilayer for the A'B case. The t' values shown in table lllll 
reflect the distance between the layers, d ~ 3.72 A for 
A'B and d = 3.60 A for AB', leading to t' = 0.25 A and 
t' = 0.91 A, respectively. When different types of atoms 
sit on the top of each other, the aligned case AA', with 
two pairs of atoms superimposed has a smaller t' than 
the AB case, which is consistent with the lower energy of 
the latter. 

Final remarks. The electronic properties of the BN 
monolayer and bilayers were calculated from first princi- 
ples. Fitting a minimal tight binding to the DFT band 
structure the parameters of the empirical model were de- 
termined. The main features of the ab-initio bands are 
well described by our minimal model. The two most sta- 
ble bilayer stackings were found to be very close in en- 
ergy, suggesting that the system can occur in Nature in 
both structures. Since the BN bilayers can be obtained 
using the exfoliation method, as in the case graphene bi- 
layers, it is conceivable that this production method can 
originate the two different and most stable stacking BN 
bilayers. 



^ K. S. Novoselov, D. Jiang, F. Schedin, T. J. Khotkevich, 
S. V. Morozov, and A. K. Geim, Proc. Natl. Acad. Sci. 
U.S.A. 102, 10451 (2005). 

^ C.R. Dean, A.F. Young, I. Meric, C. Lee, L. Wang, S. Sor- 
genfrei, K. Watanabe, T. Taniguchi, P. Kim, K.L. Shepard, 
and J. Hone Nature Nanotechnology 5, 722 (2010). 

^ Wei-Qiang Han, Lijun Wu, Yimei Zliu, Kenji Watan- 
abe and Takashi Taniguchi, Appl. Phys. Lett. 93, 223103 
(2008). 

* J. C. Meyer, A. Chuvilin, G. Algara-Siller, J. Biskupek and 

U. Kaiser, Nano Lett. 9, 2683 (2009). 
® J. H. Warner, M. H. Riimmeli, A. Bachmatiuk, B. Biichner, 

ACSNano, 4, 1299 (2010). 
^ C. Lee, Q. Li, W. Kalb, X. Liu, H. Berger, R. W. Carpik 

and J. Hone, Science 328, 76 (2010). 

R. V. Gorbachev, I. Riaz, R. R. Nair, R. Jalil, L. Britnell, 

B. D. Belle, E. W. Hill, K. S. Novoselov, K. Watanabe, T. 

Taniguchi, A. K. Geim and P. Blake, To be published. 
® L. Song, L. Ci, H. Lu, P. B. Sorokin, C. Jin, J. Ni, A. G. 

Kvashnin, D. G. Kvashnin, J. Lou, B. L Yakobson and P. 

M. Ajayan, Nano Lett. 10, 3209 (2010). 
® M. Topsakal, E. Aktiirk, and S. Ciraci, Phys. Rev. B 79, 

115442 (2009). 

° G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly 



and J. van den Brink, Phys. Rev. B 76, 073103 (2007). 
^ J. Slawinska, I. Zasada and Z. Klusek, Phys. Rev. B 81, 
155433 (2010). 

^ J. Viana Gomes and N. M. R. Peres, J. Phys.: Condens. 

Matter 20, 325221 (2008). 
^ N. Marom, J. Bernstein, J. Garel, A. Tkatchenko, E. Jose- 

levich, L. Kronik and O. Hod, Phys. Rev. Lett. 105, 046801 

(2010). 

* M. J. Rayson, P. R. Briddon Rapid iterative method 
for electronic-structure eigenproblems using localised ba- 
sis functions, Comput. Phys. Commun. 178, 128 (2008) 

^ J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. 
Lett. 77, 3865 (1996). 

^ C. Hartwigsen, S. Goedecker, J. Hutter, Phys. Rev. B 58, 
3641 (1998). 

H. J. Monkhorst, J. D. Pack, Phys. Rev. B 13, 5188 (1976). 
® J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. 

Castro Neto, Phys. Rev. Lett. 99, 256802 (2007). 
® J. Slawinska, I. Zasada P. Kosiriski and Z. Klusek, 

ArXiv:1007.3238vl 

M. Terrones, J. M. Romo-Herrera, E. Cruz-Silva, F. Lopez- 
Un'as, E. Muiioz-Sandoval, J.J. Velazquez-Salazar, H. Ter- 
rones Mat. Today 10, 30 (2007) 



